Phase diagram of the 3D Axial-Next-Nearest-Neighbor Ising model 
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The three-dimensional axial-next-nearest-neighbor Ising (ANNNI) model is studied by a modified 
tensor product variational approach (TPVA). A global phase diagram is constructed with numerous 
commensurate and incommensurate magnetic structures. The devil's stairs behavior for the model 
is confirmed. The wavelength of the spin modulated phases increases to infinity at the boundary 
with the ferromagnetic phase. Widths of the commensurate phases are considerably narrower than 
those calculated by mean-field approximations. 
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I. INTRODUCTION 

Periodically modulated magnetic structures have at- 
tracted scientific interest for several decades both exper- 
imentally and theoretically. A non-trivial phase diagram 
obtained by experimental measurements in cerium anti- 
monide (CeSb) shows a variety of different commensu- 
rately ordered magnetic structures with the underlying 
lattice P, The three-dimensional (3D) = | ax- 
ial next-nearest-neighbor Ising (ANNNI) model has been 
considered as a theoretical candidate for CeSb since it 
exhibits a rich structure when it is treated by mean- 
field approximation 3]. The 3D 5" = i ANNNI model 
is another example that shows a non-trivial spin mod- 
ulated phase — the so-called devil's stairs. This model 
has been analyzed theoretically by various approaches, 
including high-temperature series expansions 0, |^ , low- 
temperature series expansions , mean- field approxima- 
tions (21 J Monte-Carlo simulations an effective-field 
approximation (H, free-fermion methods, a phenomeno- 
logical renormalization, and other methods reviewed in 
Ref. Hillll. The Monte Carlo simulations have also been 
applied to the S = ^ ANNNI model with a finite number 
of spin layers 0|. Recently, Henkel and Pleimling con- 
sidered an anisotropic scaling at the Lifshitz point using 
the Wolff cluster algorithm and critical exponents have 
been calculated [l3| . 

The purpose of this paper is to clarify the phase struc- 
ture of the 3D 5 = i ANNNI model. Our interest 
is to study the spin modulated phases at intermediate 
temperatures, particularly, the stability of commensurate 
phases. For this purpose we apply a numerical variational 
method, the tensor product variational method (TPVA), 
to the model. In Section II we introduce the 3D ANNNI 
model and briefly discuss the variational background of 
the TPVA applied to the system. We present the numer- 
ical results in Section III where we construct the global 
phase diagram of the model and analyze the spin mod- 
ulated phases. We summarize the obtained results in 
Section IV. In Appendix, a numerical self-consistent op- 
timizing process is reviewed and efficiency of the modified 
TPVA is discussed. 



II. MODEL AND NON-UNIFORM PRODUCT 
VARIATIONAL STATE 

We study the S = ^ ANNNI model on a simple cubic 
lattice with the size L x cx) x cx) along the x, y, and z 
directions, respectively. The model is described by the 
lattice Hamiltonian 
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where the subscripts i, j, and k of the Ising spin a — ±1 
refer to the x, y, and z coordinates, respectively. The fer- 
romagnetic interaction Ji > acts between the nearest- 
neighbors and J2 > is the competing antiferromagnetic 
interaction between the next-nearest-neighbors imposed 
only in the x direction. 

Figure n shows the layer-to-layer transfer matrix T to 
the z direction which connects two adjacent spin layers 
[a] and [a] (each of the size L x 00 in the x and y di- 
rections). The transfer matrix can be exactly expressed 
as the product of partially overlapped local Boltzmann 
weights (cf. Fig.P) 
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We simplify the notations using a group of 6 spins 

{cr} = (fjjj o-j/j o-j/'j fjjj' CTj/j' o-i",j') , (3) 



with the index rule z' = z + 1, i" = z + 2, and j' = j + 1. 
The local Boltzmann weight Wf'^ of the Hamiltonian in 
Eq. ^ has the following form 
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FIG. 1: The layer-to-layer transfer matrix T[(t|(t] (left) illus- 
trated in the case for L = 5 and the local Boltzmann weight 
W^,{a\a} (right). 




FIG. 2: Graphical representation of the local variational 
weight Vij{a} (left) used to construct the trial function 
(right) in the particular case for L — 5. 



with /cb being the Boltzmann constant and the tempera- 
ture T. For reasons of simplicity and brevity, we consider 
Ji — 1 and kB—l throughout all calculations [l^ . 

We consider a variational problem for the transfer ma- 
trix T[(t|ct]. For a given trial state j^I'), the variational 
partition function per layer is given by 
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The TPVA is a numerical variational method that as- 
sumes a trial function written by the product of local 
weights V. For the ANNNI model, ^ is written in 
the product form of mutually overlapped local weights 
(cf. Fig.© 
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where we have used the simplified notation in Eq. 

In order to study non-uniform spin modulated phases, 
the local variational weights Vij{a} must be position de- 
pendent along the x direction. Each V thus contains 
2^ = 64 adjustable parameters. Since we have written 
the trial function 5" as well as the transfer matrix T in 
the product forms, both the numerator of Eq. ||SJ) 

L-2 +00 
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(7) 
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FIG. 3: The global phase diagram of the 3D ANNNI model 
obtained by the TPVA. The Lifshitz point Pl is denoted by 
the black circle. The dotted lines enclose extremely narrow 
commensurate phases. 



and its denominator 
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(*!*)= E n n (^^.w)' (8) 

[ct] 4=1 j = -00 

have also the product forms. Thus these quantities can be 
accurately calculated by means of renormalization tech- 
niques, particularly, we used the density matrix renor- 
malization group (DMRG) [ll[Tl|. 



III. RESULTS 

Figure|31shows the global phase diagram of the ANNNI 
model obtained by the TPVA. It consists of 

(i) a paramagnetic (disordered) phase, 

(ii) a uniformly ordered ferromagnetic phase, 

(iii) an antiphase with the periodic spin alignment 
(■ • • ttii ■ ■ ■) for which we use the notation (2) 
in the following, and 

(iv) a rich area of spin modulated phases. 

Region of the spin modulated phases separates the an- 
tiphase from the paramagnetic phase. The paramag- 
netic, ferromagnetic, and the modulated phases meet 
at the Lifshitz point Pl- In our calculations, it is 
located at and /cbJl/ Ji=3.83 and IS m 

agreement with the latest Monte Carlo calculations car- 
ried out by Pleimhng and Henkel J;^/Ji=0. 270(4) and 
A:bTl/Ji=3. 7475(50) (from Ref. (IJ]). 
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FIG. 4: The five selected lines AiBi, A2B2, A3B3, Pl-Bq, and F^G. 5: The spontaneous magnetization {a,) versus lattice 
A4B4 in the region of modulated phases. The points Ci and size in the x direction calculated at Ci (the upper graph) and 
C2 are marked by the white triangles at C2 (the lower one) with the lattice size L = 401. In order 

to show the data in detail, we plot i = 120, . . . , 170. 



The resulting phase diagram does not contradict to 
previous knowledge of the model. The phase bound- 
ary lines separating the ferromagnetic phase, the para- 
magnetic phase, the antiphase, and the spin modulated 
phases coincide with those obtained using the Monte 
Carlo calculations 0. We found new features of the 
model in the region of the modulated phases, where the 
Monte Carlo simulations have not yielded a satisfactory 
answer. Our results are thought of as a supplement to 
the achievements computed by the mean-field approxi- 
mations 0, at higher temperatures and by the low- 
temperature series expansions valid at low temperatures 

In the rest of this section, we focus on the region of the 
modulated phases that contains a multitude of various 
commensurate and incommensurate phases. For exam- 
ple, in Fig. 121 we plotted a few narrow areas of typical 
commensurate phases such as (3,2) = (• • ■ TTTii 

(3) = (... TTTiii •••), (4) = (••• TTTTiiii •■•), etc. 

all enclosed by the dotted lines. Note that the widths 
of these phases are substantially narrower compared to 
the mean-field approximation [3 and the effective-field 
approximation 



listed in Table I. When we obtain the spontaneous mag- 
netization, we compute the corresponding wavelengths 
by means of the Fourier transform. 

Figure El shows the spin polarizations (ct;) at the two 
parameter points: Ci on the line AiBi and C2 on A2i?2- 
These two points are chosen near the phase boundaries. 
The spin polarization at Ci exhibits the commensurate 
phase (35,2) = (••• iiiTTTiiiTTniiTT •••) with the 
wavelength A = 17/3 (the upper graph). The lower graph 
shows the commensurate phase at C2 with A « 16.7 on 
the same region along the x direction. 

Now, we give a brief discussion on influence of bound- 
ary conditions imposed to the system on the resulting 
spin polarization. In Fig. we plot (cTi) for three dif- 
ferent types of the boundary conditions. We consider a 
lattice with the size 401 x 00 x 00 and analyze the data at 
C2. On the upper graph, the spin polarization is calcu- 
lated for the fixed boundaries on the left end (the spins 
are aligned to the 'up' direction) and the free boundaries 
on the right end. The Fourier transform applied to the 
whole region i = 0, 1, . . . , 400 yields A = 16.7 ± 0.53. On 
the intermediate graph, the parallel fixed boundary con- 



A. Wavelength analysis 

We first explain relation between the conventional no- 
tation and and modulation wave length A. The antiphase 
(2) = TTiJ-TTii has periodicity of 4 lattice sites, thus 
A = 4. Another examples is the high-order commensu- 
rate phase (3, (3, 2)^) which represents the periodic spin 
sequence (TTniiTniiTT) and yields A = 26/5. 

We calculate the spin modulations along five represen- 
tative lines as depicted in Fig. ^ with their ending points 
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FIG. 6: The spontaneous magnetization obtained at C2 for 
the lattice size 401 x oo x 00. The upper, middle, and lower 
graphs display (ai) for the three different boundary condi- 
tions. 
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FIG. 7: The divergence of A the boundary with the ferromag- 
netic phase. The black circles, squares, and triangles, respec- 
tively, correspond to the selected points on the lines AiBi, 
A2B2, and A3B3. The inset shows behavior of A around the 
commensurate phase (3). The inset corresponds to the mag- 
nified area denoted by the dotted rectangle. 



ditions on both sides are imposed (the spins are aligned 
'up' at the ends). It gives A = 16.7 ± 0.41. Finally, 
the lower graph shows the anti-parallel fixed boundaries 
('up' on the left end and 'down' on the right end) with 
A = 16.7 ± 0.41. The choice of the boundary conditions 
does not affect the numerical results significantly. The 
larger lattice size is considered, the less influence of the 
boundaries is obtained, especially, off of the phase bound- 
aries. 

In Fig. Owe plot the wavelength with respect to J2/J1 
calculated on the three lines AiBi, A2B2, and A3B3. 
The dotted line is a guide for the eye to point out this 
structure. Near the boundary with the ferromagnetic 
phase, the wavelength rapidly increases. This is contra- 
dictory to the known results coming from the mean-field 
approximation. The mean-field approximation yields the 
first-order transitions between ferromagnetic phase and 
the individual commensurate phases on the boundary 
line (in details, see Ref. 0)- 

In the inset of Fig. [7| we plot details of the wave- 
length A in the vicinity of the commensurate phase (3). 
We observed that the commensurate phases (3), (3^,2), 
and (3, 2) 'lock-in' at small regions of J2/J1 (on the line 
^3,63) and the so-called "devil's stairs" behavior is ob- 
served contrary, the stairs-like structure is 
not visible at higher temperatures, as seen on the line 
AiBi near the paramagnetic boundary. 

Figure O shows A on the line from the point Bq to the 
Lifshitz point Pl- The wavelength diverges toward the 
calculated Lifshitz point at J2/J1 = 0.26. The stairs-like 
structure reveals if we zoom in on the phase diagram. For 
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FIG. 8: The behavior of the wavelength A on the line P^Bq. 
Several commensurate phases phases with the integer value 
of A are labeled. The inset shows the area around the the 
commensurate phase (5). 



this reason, we selected an area depicted by the dotted 
rectangle therein. The inset shows the magnified area 
with the commensurate phases (5) and (5,4) where the 
wavelength locks-in. 
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FIG. 9: The calculation of A on the line A^B^ at the temper- 
ature ksT/Ji = 1.5. The dotted rectangle borders an area 
shown in Fig. 1101 The inset illustrates the behavior of the 
corresponding wave vector q. 



B. Low temperature behavior 

To compare our results with analytical predictions, 
particularly, with the low-temperature series expansions 
(LTSE) 0, we have selected the line which corre- 

sponds to the temperature k^T/ Ji = 1.5. The computed 
data are plotted in Fig. |51 The commensurate phase (3) 
(A — 6) locks-in and forms a well-visible plateau. Note 
that the mean-field calculations ^J] do not exhibit any 
phases with A > 6 at k-QT/ Ji % 2. It should be also 
noted that both the mean-field approximation (at higher 
temperatures) and the LTSE (at very low temperatures) 
do not result the spring of phases with A > 6 from the 
multi-phase point J1I J\ = 0.5. 

Figure ^| illustrates the stairs- like structure of A and 
corresponds to the magnified area shown by the dotted 
rectangle in Fig. |^ The commensurate phases lock-in 
at rational values of A and are separated by the high- 
order commensurate phases and (possibly) the truly in- 
commensurate ones. Several commensurate phases are 
denoted above the stairs-like curve in Fig. ^| 

The LTSE yields a spring of an infinity of the commen- 
surate phases, such as (3) and (3, 2") for n = 1, 2, 3, . . ., 
which separate the ferromagnetic phase and the an- 
tiphase. The transition from the region of the spin mod- 
ulated phases (A > 4) to the antiphase (A = 4) does not 
contradict to LTSE. The existence of additional interme- 
diate phases, (3, 2", 3, 2"+^) n = 1,2,3,..., at a higher 
temperatures was later reported by the LTSE. Our re- 
sults contain all these commensurate phases. Moreover, 
we found unpredicted phases. For example, the tran- 
sition between the phases (3) and (3, 2) is not of the 
first order as reported by LTSE. We found that these 
two phases are separated by many commensurate phases. 



FIG. 10: The resulting devil's stairs with A > 6 observed on 
the line A4B4. 



e. g., (3", 2) and (3, (3, 2)"+i) for n = 1, 2, 3, . . . and the 
others of higher-orders as reported in Ref. We ob- 
tained such rich spin modulated structure also for A > 6, 
see Fig. Cni 

In Fig. ^] we depict our numerical results obtained at 
fceT/Ji = 1.5 and compare them with the results ob- 
tained by LTSE. The notation (7^2) corresponds to the 
antiphase wave vector q — tt/2. Note that while the 
LTSE gives the first-order transition among individual 
commensurate phases, our calculations yield subsequent 
stairs-like structures among them. Moreover, the LTSE 
calculations do not yield the commensurate phases with 
A > 6. 

We, therefore, conjecture that the 'complete' devil's 
stairs structure exists at intermediate temperatures. The 
complete devil's stairs structure suggests there are no 
first-order transition [Tol |. 

Here, we summarize those commensurate phases which 
were obtained by the numerical analysis of this model. 
Between two main commensurate phases {p) and (p+l), 
where p = 2,3,4,..., new high-order commensurate 
phases are present, such as {p"~^,p -\- 1), (p, {p + 1)"), 
{p,{p,p + 1)""^), and {{p,p + l)",j3 + 1), with n = 
2,3,4, .. . Subsequently, the following higher-order com- 
mensurate phases were found ((p"^^,p + l)"^,p",p + 1) 
and ({p, (p,p+ l)"+i)",p, (p,p+ 1)"), for m = 2, 3,4, . . . 
etc. 



IV. SUMMARY 

We applied the modified TPVA to the 3D ANNNI 
model and obtained the global phase diagram. The lo- 
cation of the Lifshitz point agrees with calculations per- 
formed by the high-temperature series expansions 0, Q 
and the recent Monte Carlo calculations The mod- 
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FIG. 11: Comparison of the numerical results at ksT/Ji = 
1.5 (the black squares) with the low-temperature series expan- 
sions represented by the dashed stairs-like curve (in Ref. '^). 



ulated phase exhibits very complex structures. We found 
that (1) the commensurate phases are substantially nar- 
rower than those reported so far, (2) the wavelength 
of the spin modulated (commensurate) phases diverges 
at the boundary with the ferromagnetic phase, (3) the 
commensurate phases merge at low temperatures tend- 
ing toward the multi-phase point J2/ Ji = 0.5 and at low 
temperatures, the wavelengths with A > 6 are obtained, 
and (4) many (possibly infinity) phases have been found 
within the modulated phase, which have not yet been 
reported. 
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Appendix 

Optimizing process for the local variational weights 

We briefly describe the optimizing process of finding 
out the local variational weights in order to maximize 
Eq. Q. This optimization is based on the self-consistent 
equation in the TPVA to achieve the minimum of the 
free energy. Numerical details in the TPVA has been 
reported in Refs. 0. 



In order to maximize the variational partition function 
in Eq. Q, by a proper tuning of the local variational 
weights V, we define two objects. One is the matrix ob- 
ject B that represents a punctured classical system 
defined on the 2-layer spin system which corresponds to 
the numerator (^ITI^*) of the variational partition func- 
tion. It is defined as 

B,,ow\a} = w^owi^} E n n^MW 

[cr],[5] ^5^0 

X Wl,{a\a}VkA^}. (9) 

Analogously, the vector object A is the punctured system 
defined on the 1-layer spin system, 

AoW-E n H^mW^^mM- (10) 

[a] 

The configuration sums in Eqs. @ and (|1U|I are taken 
over all the spin variables a except for the 6 ones at the 
center of the system. In particular, except for 



{a} — ((Ti^o cr.j+2,0 cr-j^i cr.j+i,i 0-^+2,1) 



(11) 



and analogously for {a} in Eq. @. The notations 
rife^^i rif^^o exclude Vifi{a} and Vi^o{a} from the prod- 
uct. Having defined these two objects, the variational 
partition function can be transformed into the expres- 
sion, 



E V,,o{a}B,,oW\a}VM^} 



(12) 



Now, consider a variation of Avar with respect to vari- 
ations of the local variational weights 



(5* 



E 



6X, 



SV^. 



(13) 



Carrying out the extremal condition, SX/SVi_j ~ 0, the 



self-consistent equation for the local variational weights 



Vij is then obtained 



T/ncw r^l 



f-1 ^^,0W} 



The improvement of V is performed as 



(14) 



(15) 



through 64-parameter local parameter adjust. The self- 
consistent relation, Eq. (|14|l . is a non-linear equation 
since Bi^o and Ai^o themselves depend on V. The con- 
vergence parameter s controls the rate at which the im- 
provement process of V is performed. 

Consequently, we compute the free energy per strip 



^now — — fcRTlnAv 



(16) 



and compare with the free energy J^oid calculated with 
the previous V. 
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TABLE 11: The critical temperature Tc for the 3D Ising model 
(J2 = 0). We calculate the relative errors e with respect to 
obtained by Monte Carlo simulations T^]. 



Numerical method 



e/[%] 



Mean-field approximation 6.000 33.0 

Kramers- Wannier approximation [2^ 4.587 f .7 

TPVA with 16 parameters [17] 4.570 1.3 

TPVA with 64 parameters 4.554 0.9 

Monte Carlo simulations 4.512 — 



Efficiency of the algorithm 

After the trial state [vl/} is optimized, we calculate the 
spontaneous magnetization at a site 



(17) 



Since the competing interactions exist only along the x 
direction, the system is translation invariant with respect 



to the y and z directions. Therefore, the spontaneous 
magnetization {(Ji j) is independent on j and we used 
((Ti) instead. 

In order to estimate the numerical accuracy of the im- 
proved TPVA, we compare the calculation of the criti- 
cal temperature in the pure Ising model, i. e., when 
J2 = 0, with other numerical methods. Table II summa- 
rizes the obtained Tc. It is obvious that the mean-field 
approximation overestimates Tc and does not give reli- 
able results near the phase boundaries. The improved 
TPVA with the 64 variational parameters results better 
Tc than the original TPVA with 16 parameters [HIIJ. 

We set up the convergence parameter \e\ = 10~^. As- 
suming any |e| < 10~^ is sufficient for the most cases. 

In Fig.^]we illustrate an example which demonstrates 
the systematic decay of the free energy during the DMRG 
sweeping process until it finally converges. After the 
DMRG infinite system method (ISM) is finished, the 
first left-right sweep (ILR) proceeds followed by the first 
right-left sweep (IRL) and so on. Each step of the fi- 
nite system method decreases the free energy until its 
minimum is reached. 
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